!/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!*  Copyleft (!c) 2004-2005 Div. Nucl. Med., Cyric, Tohoku Univ.    *
!*  Copyleft (!c) 2005-2011 Elseif Laboratory.                      *
!*  elseifkk@users.sf.net                                           *
!*                                                                  *
!*  All Rights Reversed.                                            *
!*                                                                  *
!*                                                                  *
!*  LICENSE AGREEMENT:                                              *
!*                                                                  *
!*  You agree that we don't agree to anything.                      *
!*                                                                  *
!*                                                                  *
!*  DISCLAIMER OF WARRANTIES:                                       *
!*                                                                  *
!*  This software is provided to you "AS IS," and the authors and   *
!*  contributors disclaim any warranty or liability obligations to  *
!*  you of any kind, implied or express, including but not limited  *
!*  to the warranties of merchantability, fitness for a particular  *
!*  purpose, and non-infringement.                                  *
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
!
! I/O routines for analize format file. 
!
! TTD:
!   1) error code and error handling routine.    
!
! note:
!   1) general image handling routines were moved to buty.
!   2) file I/O routines use some external procs written in asm,
!      in which procs use system call directly. So this utility
!      cannot work when complied on other systems than FreeBSD.
!   3) image matrix in a file, i.e., *.img will be converted 
!      to real*8 or real*4 right after it was read. 
!
module ana
  use buty
  implicit none
  ! general parameters
  integer*4,private,parameter::HEADER_SIZE    = 348
  integer*4,private,parameter::PATH_LEN_MAX   = 256
  ! ana data status
  integer*4,private,parameter::ADS_SWAPED        = Z'0001'
  integer*4,private,parameter::ADS_DT_CORRECTED  = Z'0002'
  integer*4,private,parameter::ADS_NZ_CORRECTED  = Z'0004'
  integer*4,private,parameter::ADS_CALIBRATED    = Z'0008'
  integer*4,private,parameter::ADS_TRUNCATED     = Z'0010'
  ! ana data flags
  integer*4,public,parameter::ADF_IGNORE_CF      = Z'0001'
  integer*4,public,parameter::ADF_TRUNC_POS      = Z'0002'
  ! error codes
  integer*4,parameter::ERR_READ_ANA   = 1
  integer*4,parameter::ERR_WRITE_ANA  = 2
  integer*4,parameter::ERR_READ_PGMS  = 3
  integer*4,parameter::ERR_INVALID_DATATYPE = 1
  integer*4,parameter::ERR_HEADER_SIZE      = 2
  integer*4,parameter::ERR_IMAGE_SIZE       = 3
  integer*4,parameter::ERR_ALLOCATE         = 4
  integer*4,parameter::ERR_NT               = 5
  integer*4,parameter::ERR_HEADER_ACCESS    = 6
  integer*4,parameter::ERR_IMAGE_ACCESS     = 7
  integer*4,parameter::ERR_HEADER_PROCESS   = 8
  integer*4,parameter::ERR_IMAGE_PROCESS    = 9
  integer*4,parameter::ERR_FILE_EMPTY       = 10
  integer*4,parameter::ERR_FILE_INCONSISTENT= 11
  ! these are defined in futyA.asm
  integer*4,parameter::ERRSYS_NTSTR = 1024
  integer*4,parameter::ERRSYS_OPEN  = 5
  integer*4,parameter::ERRSYS_READ  = 3
  integer*4,parameter::ERRSYS_WRITE = 4
  type ana_data
     character(PATH_LEN_MAX) fname
     integer*4 stat
     integer*1 hdr(HEADER_SIZE)
     real*8 r ! reduction factor
     real*8 c ! calibration factor (optional)
     type(image_data) id
     type(ana_data),pointer::ad
  end type ana_data
contains
  !------------------------------------------------------!
  character*32 function res_err_code(rc)
    integer*4,intent(in)::rc
    select case(rc)
    case(ERR_INVALID_DATATYPE)
       res_err_code="invalid data type"
    case(ERR_HEADER_SIZE)
       res_err_code="incorrect header size"
    case(ERR_IMAGE_SIZE)
       res_err_code="incorrect image size"
    case(ERR_ALLOCATE)
       res_err_code="allocate"
    case(ERR_NT)
       res_err_code="not implimented"
    case(ERR_HEADER_ACCESS)
       res_err_code="accessing header"
    case(ERR_IMAGE_ACCESS)
       res_err_code="accessing image"
    case(ERR_HEADER_PROCESS)
       res_err_code="processing header"
    case(ERR_IMAGE_PROCESS)
       res_err_code="processing image"
    case default
       res_err_code="not used"
    end select
  end function res_err_code
  !------------------------------------------------------!
  character*32 function res_err_system(rc)
    integer*4,intent(in)::rc
    select case(rc)
    case(ERRSYS_OPEN)
       res_err_system="open"
    case(ERRSYS_NTSTR)
       res_err_system="during null-terminate"
    case(ERRSYS_READ)
       res_err_system="read"
    case(ERRSYS_WRITE)
       res_err_system="write"
    case default
       res_err_system="undefined"
    end select
  end function res_err_system
  !------------------------------------------------------!
  character*32 function res_err_orig(rc)
    integer*4,intent(in)::rc
    select case(rc)
    case(ERR_WRITE_ANA)
       res_err_orig="write_ana"
    case(ERR_READ_ANA) 
       res_err_orig="read_ana"
    case default
       res_err_orig="unknown"
    end select
  end function res_err_orig
  !------------------------------------------------------!
  subroutine err_ana(rc)
    use cuty
    integer*4,intent(in)::rc
    integer*4 sys_rc,err_orig,err_code
    character*256 mes
    integer*4 l
    sys_rc=ishft(rc,-16)
    err_orig=ishft(iand(rc,Z'FF00'),-8)
    err_code=iand(rc,Z'FF')
    l=conc(mes,1,"Error at ")
    l=conc(mes,l,trim(res_err_orig(err_orig)))  
    l=conc(mes,l,": code=")
    l=conc(mes,l,err_code)
    l=conc(mes,l," (")
    l=conc(mes,l,trim(res_err_code(err_code)))
    l=conc(mes,l,")")
    call messn(mes(1:l-1))
    if(sys_rc/=0) then
       l=conc(mes,1,"By system call '")
       l=conc(mes,l,trim(res_err_system(sys_rc)))
       l=conc(mes,l,"'")
       call messn(mes(1:l-1))
    endif
  end subroutine err_ana
  !------------------------------------------------------!
  subroutine init_ana_data(x)
    type(ana_data),intent(out)::x
    call init_hdr(x%hdr)
    x%fname=""
    x%stat=0
    x%id%stat=0
    x%r=0
    x%c=0
    call init_image_data(x%id)
  end subroutine init_ana_data
  !
  subroutine init_hdr(x)
    integer*1,intent(inout)::x(:)
    x=0
    x( 1)=Z"5C"
    x( 2)=Z"01"
    x( 5)=Z"64"
    x( 6)=Z"73"
    x( 7)=Z"72"
    x(15)=Z"6D"
    x(16)=Z"30"
    x(39)=Z"72"
    x(40)=Z"30"
    x(41)=Z"04"
    x(49)=Z"01"
    x(57)=Z"6D"
    x(58)=Z"6D"
    x(79)=Z"80"
    x(80)=Z"40"
  end subroutine init_hdr
  !------------------------------------------------------!
  subroutine cp_ana_data(a,b)
    type(ana_data),intent(in)::a
    type(ana_data),intent(out)::b
    b%fname=a%fname
    b%stat=a%stat
    b%hdr=a%hdr
    b%r=a%r
    b%c=a%c
  end subroutine cp_ana_data
  !------------------------------------------------------!
  subroutine print_ana_hdr(hdr)
    use cuty
    integer*1,intent(in)::hdr(:)
    character*128 str
    integer*4 l
    l=conc(str,conc(str,1,"Data type: "),hdr(70+1),KIND_INTEGER2) 
    call messn(str(1:l-1))
    l=conc(str,conc(str,1,"Bpp: "),hdr(72+1),KIND_INTEGER2)
    call messn(str(1:l-1))
    l=conc(str,conc(str,1,"CF.: "),hdr(112+1),KIND_REAL4)
    call messn(str(1:l-1))
    l=conc(str,conc(str,conc(str,conc(str,1,&
         "GL. min/max: "),int(hdr(138+1),KIND=4)),"/"),&
         int(hdr(142+1),KIND=4))
    call messn(str(1:l-1))
  end subroutine print_ana_hdr
  !------------------------------------------------------!
  subroutine print_ana_data(x,title)
    use cuty
    type(ana_data),intent(in)::x
    character*(*),intent(in),optional::title
    character*128 str
    integer*4 l
    call messn("File: "//trim(x%fname))
    l=conc(str,1,"Stat:")
    if(x%stat==0) then
       l=conc(str,l," ok")
    else
       if(iand(x%stat,ADS_SWAPED)/=0) l=conc(str,l," SWAPED")
       if(iand(x%stat,ADS_DT_CORRECTED)/=0) l=conc(str,l," BAD_DT")
       if(iand(x%stat,ADS_NZ_CORRECTED)/=0) l=conc(str,l," BAD_NZ")
       if(iand(x%stat,ADS_CALIBRATED)/=0) l=conc(str,l," CAL.")
       if(iand(x%stat,ADS_TRUNCATED)/=0) l=conc(str,l," TRUNC.")
       l=conc(str,conc(str,conc(str,l," ("),x%stat,1),")")
    endif
    call messn(str(1:l-1))
    l=conc(str,conc(str,1,"RF.: "),x%r)
    call messn(str(1:l-1))
    call print_ana_hdr(x%hdr)
    call print_image_data(x%id,title)
  end subroutine print_ana_data
  !------------------------------------------------------!
  subroutine write_ana(out,output_datatype,rc)
    use futy
    use muty
    type(ana_data),intent(inout)::out
    integer*4,intent(in)::output_datatype
    integer*4,intent(out)::rc
    integer*4 s,nxyz
    integer*1,allocatable::buf(:)
    character(PATH_LEN_MAX) fname
    integer*4 istat,rc_sys
    integer*2 obp,odt
    real*8 calibration_factor
    real*8 norm
    real*8 pmax,rmin,rmax
    odt=output_datatype
    select case(odt)
    case(DT_IKIND_1)
       obp=8
    case(DT_IKIND_2)
       obp=16
    case(DT_RKIND_4,DT_IKIND_4)
       obp=32
    case(DT_RKIND_8)
       obp=64
    case(0)
       call cpmemw(odt,out%hdr(70+1),1)
       call cpmemw(obp,out%hdr(72+1),1)
       if(iand(out%stat,ADS_DT_CORRECTED)/=0) odt=get_datatype(int(obp))
    case default
       rc=ior(ishft(ERR_WRITE_ANA,8), ERR_INVALID_DATATYPE)
       return
    end select
    call cpmemw(out%hdr(70+1),odt,1)
    call cpmemw(out%hdr(72+1),obp,1)
    call real4_to_real8(1,out%hdr(112+1),calibration_factor,0.d0)
    if(calibration_factor==0) calibration_factor=1
    if(iand(out%stat,ADS_CALIBRATED)/=0) then
       norm=1/calibration_factor
    else
       norm=1
    endif
    if(out%r/=0) norm=norm/out%r ! reduction factor
    if(out%c/=0) calibration_factor=calibration_factor/out%c
    !
    call cpmemd(out%hdr(112+1),sngl(calibration_factor),1)
    call cpmemd(out%hdr(80+1),sngl(out%id%d(1)),1)
    call cpmemd(out%hdr(84+1),sngl(out%id%d(2)),1)
    call cpmemd(out%hdr(88+1),sngl(out%id%d(3)),1)
    call cpmemw(out%hdr(42+1),out%id%n(1),1)
    call cpmemw(out%hdr(44+1),out%id%n(2),1)
    call cpmemw(out%hdr(46+1),out%id%n(3),1)
    !
    nxyz=out%id%n(1)*out%id%n(2)*out%id%n(3)
    select case(out%id%datatype)
    case(DT_IKIND_1)
       rmin=minval(out%id%img_i1)
       rmax=maxval(out%id%img_i1)
    case(DT_IKIND_2)
       rmin=minval(out%id%img_i2)
       rmax=maxval(out%id%img_i2)
    case(DT_IKIND_4)
       rmin=minval(out%id%img_i4)
       rmax=maxval(out%id%img_i4)
    case(DT_RKIND_4)
       rmin=minval(out%id%img_r4)
       rmax=maxval(out%id%img_r4)
    case(DT_RKIND_8)
       rmin=minval(out%id%img_r8)
       rmax=maxval(out%id%img_r8)
    end select
    pmax=get_pmax(int(odt))
    rmax=max(abs(rmin),rmax)
    if(rmax*norm>pmax) then
       norm=pmax/rmax
       calibration_factor=1/norm
    end if
    call real8_to_real4(1,calibration_factor,out%hdr(112+1),0.d0)
    ! write header
    fname=trim(out%fname)//".hdr"
    rc_sys=fopenwrite(trim(fname)//char(0),out%hdr,HEADER_SIZE,0)
    if(rc_sys/=0) then
       rc=ior(ishft(abs(rc_sys),16),&
            ior(ishft(ERR_WRITE_ANA,8), ERR_HEADER_PROCESS))
       return
    endif
    !
    s=nxyz*obp/8
    allocate(buf(s),stat=istat)
    if(istat/=0)then
       rc=ior(ishft(ERR_WRITE_ANA,8), ERR_ALLOCATE)
       return
    endif
    !
    if(norm==1) norm=0
    select case(out%id%datatype)
    case(DT_RKIND_8)
       select case(odt)
       case(DT_IKIND_1)
          call real8_to_int1(nxyz,out%id%img_r8,buf,norm)
       case(DT_IKIND_2)
          call real8_to_int2(nxyz,out%id%img_r8,buf,norm)
       case(DT_IKIND_4)
          call real8_to_int4(nxyz,out%id%img_r8,buf,norm)
       case(DT_RKIND_4)
          call real8_to_real4(nxyz,out%id%img_r8,buf,norm)
       case(DT_RKIND_8)
          call cpmemd(buf,out%id%img_r8,nxyz*8/4)
          if(norm/=0) out%id%img_r8=nint(dble(out%id%img_r8*norm))
       end select

    case(DT_RKIND_4)
       select case(odt)
       case(DT_IKIND_1)
          call real4_to_int1(nxyz,out%id%img_r4,buf,norm)
       case(DT_IKIND_2)
          call real4_to_int2(nxyz,out%id%img_r4,buf,norm)
       case(DT_IKIND_4)
          call real4_to_int4(nxyz,out%id%img_r4,buf,norm)
       case(DT_RKIND_4)
          call cpmemd(buf,out%id%img_r4,nxyz)
          if(norm/=0) out%id%img_r4=nint(dble(out%id%img_r4*norm))
       case(DT_RKIND_8)
          call real4_to_real8(nxyz,out%id%img_r4,buf,norm)
       end select

    case(DT_IKIND_4)
       select case(odt)
       case(DT_IKIND_1)
          call int4_to_int1(nxyz,out%id%img_i4,buf,norm)
       case(DT_IKIND_2)
          call int4_to_int2(nxyz,out%id%img_i4,buf,norm)
       case(DT_IKIND_4)
          call cpmemd(buf,out%id%img_i4,nxyz)
          if(norm/=0) out%id%img_i4=nint(dble(out%id%img_i4*norm))
       case(DT_RKIND_4)
          call int4_to_real4(nxyz,out%id%img_i4,buf,norm)
       case(DT_RKIND_8)
          call int4_to_real8(nxyz,out%id%img_i4,buf,norm)
       end select

    case(DT_IKIND_2)
       select case(odt)
       case(DT_IKIND_1)
          call int2_to_int1(nxyz,out%id%img_i2,buf,norm)
       case(DT_IKIND_2)
          call cpmemw(buf,out%id%img_i2,nxyz)
          if(norm/=0) out%id%img_i2=nint(dble(out%id%img_i2*norm))
       case(DT_IKIND_4)
          call int2_to_int4(nxyz,out%id%img_i2,buf,norm)
       case(DT_RKIND_4)
          call int2_to_real4(nxyz,out%id%img_i2,buf,norm)
       case(DT_RKIND_8)
          call int2_to_real8(nxyz,out%id%img_i2,buf,norm)
       end select

    case(DT_IKIND_1)
       select case(odt)
       case(DT_IKIND_1)
          call cpmemb(buf,out%id%img_i1,nxyz)
          if(norm/=0) out%id%img_i1=nint(dble(out%id%img_i1*norm))
       case(DT_IKIND_2)
          call int1_to_int2(nxyz,out%id%img_i1,buf,norm)
       case(DT_IKIND_4)
          call int1_to_int4(nxyz,out%id%img_i1,buf,norm)
       case(DT_RKIND_4)
          call int1_to_real4(nxyz,out%id%img_i1,buf,norm)
       case(DT_RKIND_8)
          call int1_to_real8(nxyz,out%id%img_i1,buf,norm)
       end select

    case default
       ! internal image representation as integer matrix 
       ! never implimented.
       rc=ior(ishft(ERR_WRITE_ANA,8), ERR_NT)
       deallocate(buf)
       return
    end select
    fname=trim(out%fname)//".img"
    rc_sys=fopenwrite(trim(fname)//char(0),buf,s,0)
    deallocate(buf)
    if(rc_sys/=0) then
       rc=ior(ishft(abs(rc_sys),16),&
            ior(ishft(ERR_WRITE_ANA,8), ERR_IMAGE_PROCESS))
       return
    endif
    rc=0
  end subroutine write_ana
  !------------------------------------------------------!
  subroutine read_ana(file,in,internal_datatype,rc)
    use futy
    use muty 
    character*(*),intent(in)::file
    type(ana_data),intent(out)::in
    integer*4,intent(in)::internal_datatype
    integer*4,intent(inout)::rc
    character(PATH_LEN_MAX) fname
    integer*1,allocatable::buf(:)
    integer*4 s,nxyz
    integer*4 istat,rc_sys
    integer*4 idt
    integer*2 datatype,bitpix
    real*8 calibration_factor
    real*8 nz_guess
    real*8 pmax,rmax,rmin
    external swap_anahdr
    idt=internal_datatype
    call init_ana_data(in)
    in%fname=file
    fname=trim(file)//".hdr"
    if(fexist(fname)/=0) then
       rc=ior(ishft(ERR_READ_ANA,8), ERR_HEADER_ACCESS)
       return
    endif
    s=sizeoffile(trim(fname)//char(0))
    if(s/=HEADER_SIZE) then
       if(s>0) s=0
       rc=ior(ishft(abs(s),16),ior(ishft(ERR_READ_ANA,8), ERR_HEADER_SIZE))
       return
    endif
    !
    rc_sys=fopenread(trim(fname)//char(0),in%hdr(1),s,0)
    if(rc_sys/=0) then
       rc=ior(ishft(abs(rc_sys),16),&
            ior(ishft(ERR_READ_ANA,8), ERR_HEADER_PROCESS))
       return
    endif
    !
    if(in%hdr(1)/=Z'5C'.or.in%hdr(2)/=Z'1') then
       if(in%hdr(3)==Z'1'.and.in%hdr(4)==Z'5C') then
          in%stat=ior(in%stat,ADS_SWAPED)
          call swap_anahdr(in%hdr)
       endif
    endif
    !
    call cpmemw(in%id%n(1),in%hdr(42+1),1)
    call cpmemw(in%id%n(2),in%hdr(44+1),1)
    call cpmemw(in%id%n(3),in%hdr(46+1),1)
    call cpmemw(datatype,in%hdr(70+1),1)
    call cpmemw(bitpix,in%hdr(72+1),1)
    call real4_to_real8(1,in%hdr(80+1),in%id%d(1),0.d0)
    call real4_to_real8(1,in%hdr(84+1),in%id%d(2),0.d0)
    call real4_to_real8(1,in%hdr(88+1),in%id%d(3),0.d0)
    !
    IF(DATATYPE==16.AND.BITPIX==32) DATATYPE=DT_RKIND_4

    if(iand(rc, ADF_IGNORE_CF)==0) then
       call real4_to_real8(1,in%hdr(112+1),calibration_factor,0.d0)
       if(calibration_factor==0) then
          calibration_factor=1
       else
          in%stat=ior(in%stat,ADS_CALIBRATED)
       endif
    else
       calibration_factor=1
    endif

    select case(datatype)
    case(DT_IKIND_1,DT_IKIND_2,DT_IKIND_4,DT_RKIND_4,DT_RKIND_8)
       if(bitpix/=get_bitpix(int(datatype))) then
          datatype=get_datatype(int(bitpix))
          in%stat=ior(in%stat, ADS_DT_CORRECTED)
       end if
    case default
       rc=ior(ishft(ERR_READ_ANA,8), ERR_INVALID_DATATYPE)
       return
    end select
    if(idt==0) idt=datatype
    !
    fname=trim(file)//".img"
    if(fexist(fname)/=0) then
       rc=ior(ishft(ERR_READ_ANA,8), ERR_IMAGE_ACCESS)
       return
    endif
    s=sizeoffile(trim(fname)//char(0))
    nxyz=in%id%n(1)*in%id%n(2)*in%id%n(3)
    if(s/=nxyz*bitpix/8)then
       nz_guess=s/(in%id%n(1)*in%id%n(2)*bitpix/8)
       if(nz_guess-int(nz_guess)==0.and.nz_guess/=0) then
          in%stat=ior(in%stat,ADS_NZ_CORRECTED)
          in%id%n(3)=int(nz_guess)
       else
          rc=ior(ishft(ERR_READ_ANA,8), ERR_IMAGE_SIZE)
          return
       end if
    endif
    !
    allocate(buf(s),stat=istat)
    if(istat/=0) then
       rc=ior(ishft(ERR_READ_ANA,8), ERR_ALLOCATE)
       return
    endif
    !
    rc_sys=fopenread(trim(fname)//char(0),buf,s,0)
    if(rc_sys/=0) then
       rc=ior(ishft(abs(rc_sys),16),&
            ior(ishft(ERR_READ_ANA,8), ERR_IMAGE_PROCESS))
       deallocate(buf)
       return
    endif
    !
    if(iand(in%stat,ADS_SWAPED)/=0) then
       select case(bitpix/8)
       case(2)
          call swap_words(nxyz,buf)
       case(4)
          call swap_dwords(nxyz,buf)
       case(8)
          call swap_qwords(nxyz,buf)
       end select
    endif
    !
    call alloc_image_data(in%id,in%id%n,idt)
    if(iand(in%id%stat,IMG_ALLOCED)==0)then
       rc=ior(ishft(ERR_READ_ANA,8), ERR_ALLOCATE)
       return
    endif
    !
    if(trunc(int(datatype),idt))then
       select case(datatype)
!       case(DT_IKIND_1) will never happen
       case(DT_IKIND_2)
          call get_minmax_int2(nxyz,buf,rmin,rmax)
       case(DT_IKIND_4)
          call get_minmax_int4(nxyz,buf,rmin,rmax)
       case(DT_RKIND_4)
          call get_minmax_real4(nxyz,buf,rmin,rmax)
       case(DT_RKIND_8)
          call get_minmax_real8(nxyz,buf,rmin,rmax)
       end select
       pmax=get_pmax(idt)
       if(iand(rc,ADF_TRUNC_POS)/=0) rmin=rmax
       rmax=max(abs(rmin),rmax)*calibration_factor
       !
       if(rmax>pmax) then
          in%r=pmax/rmax
          calibration_factor=calibration_factor*in%r
          if(calibration_factor/=1) in%stat=ior(in%id%stat,ADS_TRUNCATED)
       end if
    endif
    !
    if(calibration_factor==1) then
       calibration_factor=0
       in%stat=iand(in%stat,NOT(ADS_CALIBRATED))
    endif
    select case(idt)
    case(DT_RKIND_8)
       select case(datatype)
       case(DT_IKIND_1)
          call int1_to_real8(nxyz,buf,in%id%img_r8,calibration_factor)
       case(DT_IKIND_2)
          call int2_to_real8(nxyz,buf,in%id%img_r8,calibration_factor)
       case(DT_IKIND_4)
          call int4_to_real8(nxyz,buf,in%id%img_r8,calibration_factor)
       case(DT_RKIND_4)
          call real4_to_real8(nxyz,buf,in%id%img_r8,calibration_factor)
       case(DT_RKIND_8)
          call cpmemd(in%id%img_r8,buf,nxyz*2)
          if(calibration_factor/=0) &
               in%id%img_r4=nint(dble(in%id%img_r4*calibration_factor))
       end select

    case(DT_RKIND_4)
       select case(datatype)
       case(DT_IKIND_1)
          call int1_to_real4(nxyz,buf,in%id%img_r4,calibration_factor)
       case(DT_IKIND_2)
          call int2_to_real4(nxyz,buf,in%id%img_r4,calibration_factor)
       case(DT_IKIND_4)
          call int4_to_real4(nxyz,buf,in%id%img_r4,calibration_factor)
       case(DT_RKIND_4)
          call cpmemd(in%id%img_r4,buf,nxyz)
          if(calibration_factor/=0) &
               in%id%img_r4=nint(dble(in%id%img_r4*calibration_factor))
       case(DT_RKIND_8)
          call real8_to_real4(nxyz,buf,in%id%img_r4,calibration_factor)
       end select

    case(DT_IKIND_4)
       select case(datatype)
       case(DT_IKIND_1)
          call int1_to_int4(nxyz,buf,in%id%img_i4,calibration_factor)
       case(DT_IKIND_2)
          call int2_to_int4(nxyz,buf,in%id%img_i4,calibration_factor)
       case(DT_IKIND_4)
          call cpmemd(in%id%img_i4,buf,nxyz)
          if(calibration_factor/=0) &
               in%id%img_i4=nint(dble(in%id%img_i4*calibration_factor))
       case(DT_RKIND_4)
          call real4_to_int4(nxyz,buf,in%id%img_i4,calibration_factor)
       case(DT_RKIND_8)
          call real8_to_int4(nxyz,buf,in%id%img_i4,calibration_factor)
       end select

    case(DT_IKIND_2)
       select case(datatype)
       case(DT_IKIND_1)
          call int1_to_int2(nxyz,buf,in%id%img_i2,calibration_factor)
       case(DT_IKIND_2)
          call cpmemw(in%id%img_i2,buf,nxyz)
          if(calibration_factor/=0) &
               in%id%img_i2=nint(dble(in%id%img_i2*calibration_factor))
       case(DT_IKIND_4)
          call int4_to_int2(nxyz,buf,in%id%img_i2,calibration_factor)
       case(DT_RKIND_4)
          call real4_to_int2(nxyz,buf,in%id%img_i2,calibration_factor)
       case(DT_RKIND_8)
          call real8_to_int2(nxyz,buf,in%id%img_i2,calibration_factor)
       end select

    case(DT_IKIND_1)
       select case(datatype)
       case(DT_IKIND_1)
          call cpmemb(in%id%img_i1,buf,nxyz)
          if(calibration_factor/=0) &
               in%id%img_i1=nint(dble(in%id%img_i1*calibration_factor))
       case(DT_IKIND_2)
          call int2_to_int1(nxyz,buf,in%id%img_i1,calibration_factor)
       case(DT_IKIND_4)
          call int4_to_int1(nxyz,buf,in%id%img_i1,calibration_factor)
       case(DT_RKIND_4)
          call real4_to_int1(nxyz,buf,in%id%img_i1,calibration_factor)
       case(DT_RKIND_8)
          call real8_to_int1(nxyz,buf,in%id%img_i1,calibration_factor)
       end select

    case default
       ! internal image representation as integer matrix 
       ! will never implimented.
       rc=ior(ishft(ERR_READ_ANA,8), ERR_NT)
       deallocate(buf)
       return
    end select
    deallocate(buf)
    rc=0
  end subroutine read_ana
  !
  integer*4 function get_bitpix(dt)
    integer*4,intent(in)::dt
    select case(dt)
    case(DT_IKIND_1)
       get_bitpix=8
    case(DT_IKIND_2)
       get_bitpix=16
    case(DT_IKIND_4,DT_RKIND_4)
       get_bitpix=32
    case(DT_RKIND_8)
       get_bitpix=64
    case default
       get_bitpix=0
    end select
  end function get_bitpix
  !
  integer*4 function get_datatype(bpp)
    integer*4,intent(in)::bpp
    select case(bpp)
    case(8)
       get_datatype=DT_IKIND_1
    case(16)
       get_datatype=DT_IKIND_2
    case(32)
       get_datatype=DT_IKIND_4
       !       get_datatype=DT_RKIND_4 ! can not resolve !!!
    case(64)
       get_datatype=DT_RKIND_8
    case default
       get_datatype=0
    end select
  end function get_datatype
  !
  logical*4 function trunc(in,out)
    integer*4,intent(in)::in,out
    if(out>=in) then
       trunc=.false.
    else 
       trunc=.true.
    end if
  end function trunc
  !
  real*8 function get_pmax(idt)
    integer*4,intent(in)::idt
    select case(idt)
    case(DT_IKIND_1)
       get_pmax=127
    case(DT_IKIND_2)
       get_pmax=32767
    case(DT_IKIND_4)
       get_pmax=2147483648
    case(DT_RKIND_4)
       get_pmax=huge(0.)
    case(DT_RKIND_8)
       get_pmax=huge(0.d0)
    case default
    end select
  end function get_pmax
  ! 
  subroutine read_pgms(file,in,nz1,nz2,rc)
    use futy
    use muty 
    use cuty
    character*(*),intent(in)::file
    type(ana_data),intent(out)::in
    integer*4,intent(in)::nz1,nz2
    integer*4,intent(inout)::rc
    character(PATH_LEN_MAX) fname
    integer*4 nx,ny,ix,iy,iz,nz
    integer*4 rc_sys
    character*32 str
    integer*4 s,s0
    integer*4 k,l,istat
    integer*4,parameter::PGM_HEADER_SIZE=74
    integer*1 hdr(74)
    integer*1,allocatable::img(:)
    integer*4 nx1,nx2
    logical*4 interpolate,first
    rc=0
    s0=0
    nx1=0
    nx2=0
    call init_ana_data(in)
    nz=nz2-nz1+1
    interpolate=.false. ! in case some slices are missing 
    first=.true.
    do iz=nz1,nz2
       l=d2str(iz,str)
       fname=trim(file)//str(1:l)//".pgm"
       s=sizeoffile(trim(fname)//char(0))
       rc_sys=fopenread(trim(fname)//char(0),hdr(1),74,0)
       if(rc_sys/=0.or.s==0) then
          if(interpolate) then
             rc=ior(ishft(abs(rc_sys),16),&
                  ior(ishft(ERR_READ_PGMS,8), ERR_HEADER_PROCESS))
             return
          else
             interpolate=.true.
             cycle
          end if
       endif
       if(s0==0)then
          s0=s
          str(1:1)=char(hdr(63))
          str(2:2)=char(hdr(64))
          str(3:3)=char(hdr(65))
          str(4:4)=char(hdr(66))
          read(str(1:4),*,iostat=istat) nx
          !
          str(1:1)=char(hdr(67))
          str(2:2)=char(hdr(68))
          str(3:3)=char(hdr(69))
          str(4:4)=char(hdr(70))
          read(str(1:4),*,iostat=istat) ny
          if(nx>ny)then          
             nx1=(nx-ny)/2
             nx2=nx-(ny+nx1)
             nx=ny
          endif
       else
          if(s==0) then
             rc=ior(ishft(abs(s),16),&
                  ior(ishft(ERR_READ_PGMS,8), ERR_FILE_EMPTY))
             return
          else if(s/=s0) then
             rc=ior(ishft(abs(s),16),&
                  ior(ishft(ERR_READ_PGMS,8), ERR_FILE_INCONSISTENT))
             return
          endif
       end if
       !
       if(.not.allocated(img)) allocate(img(s))
       rc_sys=fopenread(trim(fname)//char(0),img(1),s,0)
       if(rc_sys/=0) then
          rc=ior(ishft(abs(rc_sys),16),&
               ior(ishft(ERR_READ_PGMS,8), ERR_IMAGE_PROCESS))
          return
       endif
       !
       if(first)then
          in%id%n(1)=nx
          in%id%n(2)=ny
          in%id%n(3)=nz
          call alloc_image_data(in%id,in%id%n,DT_IKIND_2)
          first=.false.
       end if
       !
       k=PGM_HEADER_SIZE+nx1
       do iy=1,ny
          do ix=1,nx
             k=k+1
             if(img(k)>=0) then
                in%id%img_i2(ix,ny-iy+1,nz-iz+1)=img(k)
             else
                in%id%img_i2(ix,ny-iy+1,nz-iz+1)=256+img(k)
             endif
          end do
          k=k+nx2+nx1
       end do
       !
       if(interpolate) then
          do iy=1,ny
             do ix=1,nx
                in%id%img_i2(ix,ny-iy+1,(nz-iz+1+1))=&
                     (in%id%img_i2(ix,ny-iy+1,(nz-iz+2+1))&
                     +in%id%img_i2(ix,ny-iy+1,nz-iz+1))/2
             end do
          end do
          interpolate=.false.
       end if
    end do
    deallocate(img)
    rc=0
  end subroutine read_pgms
  !
  subroutine read_shimadzu(file,in,rc)
    use futy
    use muty 
    use cuty
    character*(*),intent(in)::file
    type(ana_data),intent(out)::in
    integer*4,intent(inout)::rc
    integer*4,parameter::SMZ_HEADER_SIZE=12288
    integer*1 hdr(SMZ_HEADER_SIZE)
    integer*4 rc_int
    integer*4 nt,it,offset
    integer*2,allocatable::buf(:)
    real*8,allocatable::rbuf(:,:,:)
    integer*4 npixel,image_size
    integer*4 iy,iz
    real*4,allocatable::frame_period(:)
    call init_ana_data(in)
    rc_int=fopenread(file,hdr,SMZ_HEADER_SIZE,0)
    call cpmemd(in%id%n(3),hdr(40+1),1)    
    in%id%n(1:2)=128
    call cpmemd(nt,hdr(44+1),1)    
    call real4_to_real8(1,hdr(1920+1),in%id%d(1),0.d0)
    call real4_to_real8(1,hdr(1924+1),in%id%d(2),0.d0)
    call real4_to_real8(1,hdr(1928+1),in%id%d(3),0.d0)
    offset=SMZ_HEADER_SIZE
    call alloc_image_data(in%id,in%id%n,DT_RKIND_8)
    npixel=product(in%id%n)
    image_size=npixel*2 ! data type = integer*2 
    allocate(buf(npixel))
    allocate(rbuf(in%id%n(1),in%id%n(2),in%id%n(3)))
    allocate(frame_period(nt))
    call cpmemd(frame_period,hdr(1232+1),nt)
    if(nt>1) then
       frame_period=frame_period/sum(frame_period)
    else
       frame_period=1
    end if
    in%id%img_r8=0
    do it=1,nt
       rc_int=fopenread(file,buf,image_size,offset)
       call int2_to_real8(npixel,buf,rbuf,0.d0)
       in%id%img_r8=in%id%img_r8+rbuf*frame_period(it)
       offset=offset+image_size
    end do
    do iz=1,in%id%n(3)
       do iy=1,in%id%n(2)/2
          call cparray(rbuf(1,1,1),in%id%img_r8(1,iy,iz),in%id%n(1))
          call cparray(in%id%img_r8(1,iy,iz),in%id%img_r8(1,in%id%n(2)-iy+1,iz),in%id%n(1))
          call cparray(in%id%img_r8(1,in%id%n(2)-iy+1,iz),rbuf(1,1,1),in%id%n(1))
       end do
    end do
    deallocate(frame_period)
    deallocate(rbuf)
    deallocate(buf)
    rc=0
  end subroutine read_shimadzu
end module ana
